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Abstract 

For a spatial characteristic, there exist commonly fat-tail frequency distributions of 
fragment-size and -mass of glass, areas enclosed by city roads, and pore size/volume 
in random packings. In order to give a new analytical approach for the distributions, 
we consider a simple model which constructs a fractal-like hierarchical network 
based on random divisions of rectangles. The stochastic process makes a Markov 
chain and corresponds to directional random walks with splitting into four particles. 
We derive a combinatorial analytical form and its continuous approximation for the 
distribution of rectangle areas, and numerically show a good fitting with the actual 
distribution in the averaging behavior of the divisions. 
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1 Introduction 



Historically, a discussion for the origin of skew distribution that appears in 
sociological, biological, and economic phenomena goes back to the Simon's 
stochastic model [1]. At the beginning of this century, the evidences of fat-tail 
degree distribution have been also observed in many real networks [2] through 
computer analyses for large data. Some cases look like a power-law distribution 
while some other cases a lognormal distribution, it is difficult to discriminate 
these distributions in general. Since the tail in a lognormal distribution re- 
semble power-law behavior, only the part may be observed. Moreover, the 
generation mechanisms of the distributions are not exclusive but intrinsically 
connected in preferential attachment, multiplicative process, and other mod- 
els [3]. On a stochastic process of geometric Brownian motion, a model of 
double Pareto distribution generates a lognormal body and Pareto tail in the 
continuous distribution [4] [5]. On other process of iterative division of cells, 
the frequencies experimentally follow such distributions of fragment-size and 
-mass of glass [6] [7], areas enclosed by city roads [8] [9] [10], pore size/volume in 
random packings [11] [12] , and areas enclosed by edges in models of urban street 
patterns [13] and geographical networks [14]. In addition, the distribution in 
fragmentation of glass changes from a lognormal to a power-law-like according 
to low and high impacts [6], which determine the limitation of breakable sizes. 
It is worthwhile to study a mechanism in abstract models for generating the 
similar distributions in spite of different physical quantities and operations in 
a variety of research fields: socio-economic, material, computer, and physical 
sciences, even if we apart from the reproduce of realistic phenomena and do 
not insist on the detail process or the macroscopic properties exactly in broken 
fragmentation of glass. 

For crack patterns, a random tessellation model has been proposed [15]. It is 
based on a stochastic point process, consisting of the division of a randomly 
chosen face (cell) according to its life-time by adding a random line segment. 
The distribution of tessellations is invariant for an appropriate rescaling, whose 
characteristic is called stable with respect to iteration (STIT). In addition, the 
length distribution of segments is analytically obtained in a classification of 
several types of segments [16] [17]. However, the distribution of areas enclosed 
by segments is not derived. The adjustment of life-time for each cell is also 
not easily applicable even in the sophisticated mathematical model of STIT, 
when we consider a construction of network in a procedural manner such as 
load balancing in a territory of node. 

Thus, one of the issues is a theoretical analysis for the distribution of areas in 
a fractal-like structure. The presence of hierarchy and scaling law is important 
[18] for understanding the universal mechanism to generate such a structure. 
Through a mathematical model, we focus on the spatial phenomena based on 
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dividing by four with a randomness for the simplicity with analogous struc- 
tures to road networks. However, for the iterative division processes, only a few 
mathematical models are known. In random tessellations, although the cell- 
selection and cell-division rules are classified into equally-likely, area-weighted, 
perimeter-weighted, corner-weighted, and so on, the length distributions in 
one-dimension are merely analyzed for some simple rules [19]. In a quadtree 
model characterized as a typical road network, the shortest paths and maxi- 
mum flow are analyzed [20], but the distribution of areas is not discussed. 

On the other hand, we recently proposed multi-scale quartered (MSQ) net- 
works based on a self-similar tiling by equilateral triangles or squares [21,22]. 
This model is constructed by iterative division of faces, and is also suitable for 
the analysis of depth distribution of layered areas in a framework of Markov 
chain [23]. Moreover, from an application point of view, the MSQ networks 
without hub nodes have several advantages of the strong robustness of con- 
nectivity against node removals by random failures and intentional attacks, 
the bounded short path as t = 2-spanner [24], the efficient face routing by 
using only local information, and a scalable load balancing performed by the 
divisions of the node's territory for increasing communication or transporta- 
tion requests. However, due to the self-similar tiling in the scalably growing 
MSQ networks, the position of a new node is restricted on the half-point of an 
edge of the chosen face, and the link length is proportional to {\) H where H 
is the hierarchical depth number of divisions. The restriction is unnatural for 
many division processes in physical or social phenomena. Thus, we generalize 
the divisions of squares to ones of rectangles with any link lengths instead of 
the iterative halvings. 

The organization of this paper is as follows. In Sec. 2, we introduce a gener- 
alization of MSQ network model. In Sec. 3, for the distribution of areas, we 
derive the exact solution on a combinatorial analysis. We point out that the 
behavior of divisions is equivalent to directional random walks with splitting 
into four particles. The representation of random walks with splitting gives us 
an inspiration for the combinatorial analysis. However, this approach is lim- 
ited for application to very small networks. In Sec. 4, we consider a continuous 
approximation of the distribution of areas for large networks. We decompose 
the distribution function into two components of Poisson and gamma distri- 
butions, and numerically investigate the fitting of the mixture distribution 
by the two components with the actual distribution of areas in the divisions 
for generating a fractal-like network structure. In Sec. 5, we summarize these 
results. 
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2 Division process 



We consider a two-dimensional LxL square, whose lattice points (0, 0), (0, 1), 
. . ., (0, L), . . ., (L,0), (L, 1), . . ., (L, L) in the x-y coordinates give the feasible 
setting positions of nodes. Initially, there exist only the outer square with 
four corner nodes and edges of the length L. The following model generalizes 
the MSQ network [21,22,23] from recursive divisions of squares to ones of 
rectangles. In the MSQ network based on a self-similar tiling by squares, the 
division is restricted at the half point of an edge, therefore we extend the 
division to that at a cross point of the vertical and horizontal segments on the 
square lattice. The case of L — > oo gives a general position for the division 
point. Varying the value of L controls the limitation of divisible size relatively, 
as similar to low and high impacts in the fragmentation of glass [6]. 

The proposed network is iteratively constructed for a given L as follows. At 
each time step, a rectangle is chosen uniformly at random (u.a.r), and it is 
divided into four smaller rectangles. Then, the smaller rectangle with an area 
x x y (x,y denote the two edge lengths) is generated from the chosen rectangle 
with an area x' x y' . Simultaneously, rectangles with the areas (x' — x) x y, 
xx (y' — y), and (x' — x) x (y' — y) are generated. Here, two division axes are 
chosen u.a.r from the horizontal and vertical segments of an L x L lattice (see 
the left of Fig. 1). In other words, each edge length x, y e Z + — {1, 2, . . .} is 
randomly chosen as a positive integer in x+1 < x' < L and y+1 < y' < L. The 
stochastic network generation makes a Markov chain. The state is represented 
by a vector (nn, . . . , n xy , . . . , till), where n xy denotes the number of rectangle 
with the area x x y (The Markov chain is degenerately simplified by ignoring 
the difference of areas in subsection 4.1 for discussing the distribution of layers 
defined by the depth of divisions.). The stochastic process is characterized by 
the fact that the transition probability to divide a rectangle with the area 
x x y is not fixed but proportional to n xy because of the uniformly at random 
selection of a face. In other words, the probability depends on a sequence 
of chosen rectangles during the transition until a final absorbing state for the 
indivisible width x — 1 or y — 1. We remark that the minimum edge length x = 
1 or y = 1 bounds the number of the states finitely, while the MSQ networks 
[23] have the infinite states without a limitation of the subdivision. The scaling 
relation of the maximum iteration time is numerically obtained as T max ~ 
L 1 - 91 in the averaging of absorbing states. Although this paper discusses a 
simple case of uniformly random selections of rectangles and of division points 
in order to deduce the distribution of areas, the selections can be extended 
to other cases, e.g., according to a given population in a territory of node 
for real statistical data. Note that the above process is different from the 
Galton- Watson type branching process with a time-independent probability 
for generating offspring [25]. 
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Division of Rectangle Directional Random Walk 

Fig. 1. Correspondence between the division of a rectangle and the directional ran- 
dom walks with splitting. Vertical and horizontal thin lines at even intervals are the 
segments of an L x L lattice. 

3 Combinatorial analysis 



As shown in Fig. 1, the recursive generation process can be regarded as direc- 
tional random walks of particle with splitting into four copies in the framework 
of two-dimensional cellular automaton (CA), when a pair (x, y) of edge lengths 
of a rectangle is corresponded to the position of particle in the x-y coordinates. 
A particle is randomly chosen at a time step, and moves toward smaller co- 
ordinate values from (x', y') to (x, y), where x < x' and y < y', until reaching 
the boundary at x — 1 or y — 1. We emphasize that this type of CA with 
splitting differs from the asymmetric simple exclusion process (ASEP) [26] 
and the contact (or voter) model [25], since there are no spatial exclusions 
and no interaction between any particles. 

This representation of random walks with splitting is inspirational for deriving 
a combinatorial analytic form of the distribution of areas. We consider the 
number n xy of particles at {x,y), equivalently the number of rectangles with 
the area x x y. Remember that x and y are integers. The average behavior is 
described by the following system of difference equations for 2 < x,y < L — 1 
with the sum by taking over the integers x + 1 < x' < L and y + 1 < y' < L, 

A "-» = -"-» + g (*--w-ir (1) 



where An xy is the average difference of n xy in one step, and p xy = n xy / J2 x ">iy'">i n x"y" 
is the existing probability of a particle at (x,y). The factor 4 in the numer- 
ator of right-hand side of Eq.(l) is due to feasible positions of the x x y at 
left/right and upper/lower corners in the division of x' x y' . The denominator 
is the combination number for the relative positions of emanating particles 
in the intervals [1, x' — 1] and [l,y' — 1], and is equivalent to the number for 
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selecting two axes in the division of rectangle with the area x' x y'. 



From An xy = in Eq.(l), we derive 



Pl-il-i 



(i-ir 



PxL-l = PL-ly = ^ _ X>l,y>\, 
PL-2L-2 = 1 + 



(l-2) 2 ; (L — l) 2 

In general, we obtain the solution by applying the above in decreasing order 
of x and y recursively. 



PXV ~ {' + ? ( nLl to " 4, - 1)) } (^) 2 ' 



(2) 



where X)-p denotes the sum for a set of paths through the points (xi, yi), (x2, 1/2), • • • , Hi), 
with Xi, yi G Z + , x < xi < X2 < ■ ■ ■ < Xi < . . .xi < L — 1, and y < yi < y 2 < 
. . . < yi < . . . yi < L — 1 in all combinations of / = 1, 2, . . . min{L — 1 — x, L — 

i-y}- 

By substituting the solution p xy of Eq.(2) into the following right-hand sides, 
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n>xi — J2 X ': 
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-1) 



n n — Sx'>iy>i 7^ 



we obtain the distribution P(A) of rectangles with the area A. The sum is 
taken over the positive integers x' , y' < L, riu, n xl , and n ly denote the numbers 
of the finally remaining rectangles with the areas 1 x 1, x x 1, and 1 x y, 
which are no more divisible. Note that the unknown factor pll disappears 
by the numerator and the denominator in all of -P(l) = nu/M and P(x) = 
(n xl + n lx )/M, where Af = riu + J2 x ">i n x"i + J2 y ">i n iy" denotes the total 
number of the divided rectangle faces. 

Figure 2 shows the distribution of areas with width one. Our solution denoted 
by lines is almost completely fitting with the actual distribution denoted by 
marks. The part of linear tail becomes longer, as L is larger. Note that these 
distributions in Fig. 2 are estimated well by lognormal functions in this entire 
range of A. 
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Area of rectangle: A 



Fig. 2. (Color online) Distribution of areas: lx 1,1x2,... , 1 x L — 1 in the extreme 
rectangles for L = 10, 15, 20, and 25 from left to right. The (plus, cross, star, and 
rectangle) marks show the averaged result in 100 samples of the actual divisions, 
and the corresponding (red, green, blue, and magenta) lines show the solution of 
Eq.(2) on the combinatorial analysis. The short (cyan) segment guides the slope of 
—2 corresponded to the exponent in broken fragments of glass [6] and city roads 
[8,9]. 

4 Continuous approximation 

Since it possibly causes a combinatorial explosion to calculate the extreme 
distribution of areas: 1 x 1, 1 x 2, . . . , 1 x L - 1 at the absorbing states in the 
Markov chain, the application of Eq.(2) is restricted for a very small L. In 
order to analyze the distribution of areas for a large L, we approximate the 
process to be divisible at any positions on two edges of a rectangle, by ignoring 
the restriction to the segments on an L x L square lattice. For 1 = 1,2,... until 
the maximum layer at a given time, we consider the sum of the product of pi 
and g2i()ogA), which denote the frequencies of layer / and of area A in the 
layer I. Here, the number I represents the depth of divisions. We derive these 
frequencies separately. Numerical simulations in subsection 4.3 show a good 
fitting of the mixture distribution J2i PW2i(log A) with the actual distribution 
of areas in the average behavior of the divisions. 

4-1 Distribution of layers 

In this subsection, we ignore the difference of areas in each layer, and treat 
only the number of faces. Then, the stochastic division process characterized 
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as a Markov chain is simplified. Figure 3 shows the state transitions in the 
first few steps. We consider the number ni(t) of faces in the l-th layer at a 
time step t, and derive approximative solutions for the existing frequency of 
faces in the layer. 



(3,4,0 ) 



(2,8,0,....) (3,3,4,0 ) 

2/10^/ 3/10 | 



(1,12,0,....) 



(2,7,4,0,....) (3,2,8,0,....) (3,3,3,4,0,. 



Fig. 3. Branching tree diagram of the state vector (m, ri2, • • •) for the division process 
at t = 2, 3,4 steps (from top to bottom). Here, ni is the number of faces on the Z-th 
depth, when the difference of areas is ignored in the count. Each fraction denotes 
the transition probability. Note that two transitions until t = 2 through (4,0, . . .) 
at t = 1 and the initial square are trivial. 

As shown in [23] , the averaging behavior of difference 

Ani^ni{t + l)-m{t), (3) 

can be written to 

An t = mp,_ 1 (f) -Pi(t), (4) 



since a face in the layer / chosen with the probability pi(t) is divided into 
m = 4 smaller ones which belong to the layer 1 + 1, therefore a face in the 
layer I — 1 contributes to increase the number of faces in the layer /. Note that 
the simultaneously created m faces at a time belong to a same layer even with 
different areas. For a large t, by noticing n;(t) = J\f(t)pi(t) and substituting 
M{t) = J2i n i — 1 + ( m — 1)* ~ ( m — 1)* m to the right-hand side of Eq. (3), 
it is 

A ni = (m - l)(t + l) P i(t + 1) - (m - l)t Pl (t), 

= {m- l)t[ Pl {t + 1) - Pl (t)} + (m - l)p,(f + 1). 
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Using pi(t + 1) w because oft + l^t3>l, Eq. (4) is rewritten to 

Tfl 

p,(* + 1) ~Pi(t) = - {m _ 1)t {Pi( t ) ~Pi-i(t)h (5) 



where p = is assumed for convenience. The solution of Eq.(5) is not eas- 
ily derived even from a formal representation by using a generating function 
because of the combinatorial explosion involved with very complicated recur- 
sive operations [23]. For the m = 2 divisions, the difference equation (5) is 
equivalent to Eq.(14) in [19] for a crack model, however it differs to consider 
a rescaling length factor at each time step and to analyze a cumulative distri- 
bution in the one-dimensional case. 

On the other hand, we also derive the following expression [23] by using a 
model in an interacting infinite particle system, 



^ = mn i _i-n i , I > 2, (6) 

dni , 

" (7) 

The solution is 

M(t) = e {m ^\ (9) 

(mr)' 1 mT . , 

Pl ~ JT^Tji ' (10) 



where r > has a logarithmic timescale from the relation l + (m— l)t = e^ m ^ T 
of the total number of faces. Note that the distribution of Eq.(lO) coincides 
with the solution of Eq.(5) asymptotically after a huge time [23]. This form 
of ni in Eq.(8) can be applied for calculating a fractal dimension at a proper 
time in the MSQ network based on a self-similar tiling, and extended to a 
preference model for selecting a face in Appendixes A and B. 



4-2 Distribution of areas in the l-th layer 



Remember that, at each time step, a rectangle is chosen uniformly at random 
(u.a.r). For the division, vertical and horizontal axes are also chosen u.a.r from 
the segments on an L x L square lattice. We focus on a set of rectangle faces 
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on only the Z-th layer generated in some steps. After /-time selections, the 
subdivided face belongs to the layer I. The area Si is given by the product of 
shrinking rates < X; L , Yi < 1, i — 1, 2, . . . , I, for two edges of rectangle, 

Si = H^XiYiL 2 , 



where X,- L and Yi are rational numbers, and S\ is a positive integer in the 
division process, strictly speaking. 

As an approximation for a large L, we assume that the random variables X, 
and Yi follow a (0, 1) uniform distribution. Then we define a variable x = f 
— \og(Si/L 2 ) = — X)i(logXj + logFj) in the range of equivalent relation x > 
L 2 > Si, the probability of x follows a gamma distribution 

x 2l-l 

9»W = ^ X ^Y)V (") 



Here, the mean /x and the variance a 2 of log Xj are 
i 

/i = 1 logXrfX = -1, 
o 

i 

a 2 = J (logX -fi) 2 dX = 1. 



Therefore, by a central-limit theorem, we have a normal distribution N(0, 1) 
asymptotically 



log^-(logL 2 -2l) 
o\'l 



N(0,1), 



log^-2/ = log(- 



In other words, the average shrinking rate of edge is 1/e for each division of a 
rectangle, and it is slightly smaller than 1/2 for that of a square in the MSQ 
network. We can easily transform g 2 i (x) to the function of log Si by using the 
shift of x = 2 log L — log Si for a constant L. 
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4-3 Simulation for the mixture distribution 

In numerical simulation, we investigate the distributions decomposed into the 
approximative pi and g2i, and discuss the condition for a good fitting to each of 
components in the actual distributions for the divisions of rectangles. Indeed, 
our approximation of J2iPi92i shows reasonable agreement with the actual 
distribution of areas. Here, we consider the complementary cumulative distri- 
bution (CCD) of areas with a merit of smoothing effect, because the frequency- 
distribution itself has a huge variety of areas especially for a large L, therefore 
it is practically impossible to gather these samples in a proper frequency. 

Figure 4 shows the distribution pi of layers at time steps t = 50, 500, and 
5000 < T max ~ L im for L = 10 8 and 10 5 . By the effect of width one, there 
exists a gap between the solution of difference equation(5) denoted by open 
marks and the actual distribution denoted by lines, although these distri- 
butions almost coincide in the MSQ networks based on a self-similar tiling 
[23]. The gap becomes slightly larger as L is smaller in Fig. 4(b) and t is 
larger further to the right, because the effect tends to appear in more coarse- 
grained divisions and a deeper layer. With the growing of the total number 
M = 1 + (to — l)t of faces, the depth of a face tends to be deeper as the 
time step t is larger. In addition, we note the expectation (I) ~ logt [23]. The 
Poisson distribution of Eq.(10) denoted by closed marks has a slightly larger 
gap than the solution of difference equation(5) denoted by open marks. In the 
semi-log plots of Fig. 4(c) (d), we also investigate the fitting of the tails. In the 
tails, the discrepancy between the actual distribution and our approximation 
appears for the large L = 10 8 as t is larger further to the right, while these 
distributions do not fit any longer for the small L = 10 5 . Here, in order to keep 
the accuracy of approximation, we set the initial condition of {pi(0)} by the 
existing frequency of faces in each layer at t — 5 that is explicitly determined 
from the branching tree diagram with the transition probability as shown in 
Fig. 3. Note that the calculation becomes more complex as t is larger, although 
a higher accuracy is expected. 

Figure 5 shows the CCD of x = \og(L 2 /A) restricted in the /-th layer. We 
chose the most observable layer I = 5,8, 11, and 14, which correspond to the 
peaks of p\ at t — 50,500,5000, and 50000 <C T max , respectively. Because 
the most observable layer is dominant in the mixture distribution Y^iVidn- 
The effect of width one tends to appear as t is larger as shown in more right 
curves. The actual distributions denoted by lines and the gamma distribution 
g2i{x) of Eq.(ll) denoted by marks almost coincide until increasing around 
t = 5000 step (the 3rd curve from left) in Fig. 5(a) for L = 10 8 , however 
the distributions begin to differ from larger than t = 500 step (the 2nd curve 
from left) in Fig. 5(b) for L = 10 5 . Thus, the effect of width one becomes 
stronger, as L is smaller, the discrepancy between the actual distribution and 
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Fig. 4. (Color online) Distributions pi of layers at time steps t = 50, 500, and 5000 
from left to right. The (black) solid and (green) dashed lines correspond to the 
averaged result by 100 samples of the actual divisions of rectangles for (a) L = 10 8 
and (b) L = 10 5 , respectively. The semi- log plots are shown in (c) and (d) for 
investigating the fitting in the part of tail. The open and closed marks correspond 
to the solution of Eq.(5) and Poisson distribution in Eq.(10). The discrepancy of 
positions between the open marks and the lines is due to the effect of width one in 
the extreme rectangles. Note that these marks for L = 10 8 and 10 5 are the same at 
each time step, and that only the solid and dashed lines are different in (a) and (b) 
or (c) and (d). 

the approximative g2i is unignorable. We also confirm this phenomenon for 
the distributions in a commonly existing layer as shown in Fig. 6. The small 
cases of L — 10 5 and 10 6 give inaccurate approximations even at t = 500. 
In other word, for a smaller L at more coarse-grained divisions, the effect of 
width one already appears before I = 12 in a smaller (shallower) layer. This 
result is consistent with the average shrinking rate 1/e of edge per division, 
e 12 > 10 5 and e 14 > 10 6 , as mentioned in subsection 4.2. 



In order to grasp the tendency of approximation as a whole, we investigate 
the maximum difference in the CCD of areas between the actual data in the 
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Fig. 5. (Color online) CCDs of x == log(L 2 /^4) restricted in the layer / at time steps 
t = 50, 500, 5000, and 50000 from left to right, which correspond to I = 5, 8, 11, and 
14 at the peaks of pi (see Fig. 4) in this order. The (orange, red, blue, and green) 
lines show the averaged results for 100 samples of the actual divisions of rectangles, 
while the corresponding marks (square, triangle, circle, and plus from left to right) 
show the results for the gamma distribution of Eq.(ll). 
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Fig. 6. (Color online) Comparison of the CCDs of x in a commonly existing layer 
I = 12 to the sizes of I at a time step t = 500. From left to right, the (green, 
magenta, blue, and black) dashed lines show the averaged results for 100 samples of 
the actual divisions of rectangles for L = 10 5 , 10 6 , 10 7 , and 10 8 . The furthest right 
(red) solid line shows the result for the gamma distribution of Eq.(ll). 

divisions of rectangles and our approximation for several time steps. Figure 7 
shows that the difference increases as t is larger and L is smaller. In each size of 
L, the lines with open marks give slightly better approximation than the lines 
with closed marks, where the former use the solution of difference equation (5) 
and the latter use the Poisson distribution (10) as pi. The accuracy depends 
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Fig. 7. (Color online) Maximum difference in the CCDs P(> A) between the aver- 
aged result by 100 samples of the actual divisions of rectangles and our approxima- 
tion by the mixture distribution Y^i pig2i at several time steps t. From bottom to top, 
the (black, magenta, and cyan) lines with marks show the results for L = 10 8 , 10 6 , 
and 10 4 . The open and closed marks correspond to the mixtures with the solution 
of Eq.(5) and with the Poisson distribution of Eq.(10), respectively. 



on the pi as shown in Fig. 4 at a small t, while the effect of gn is added at a 
large t as shown in Fig. 5. The saturated behavior in Fig. 7 is not important, 
because the difference between any CCDs is bounded in [0, 1] by nature. 

In more detail, we investigate the CCDs of areas at t = 50 step in Fig. 8. We 
obtain a good fitting for L = 10 8 in Fig. 8(b) (d), but remark on a small gap 
for L = 10 3 in Fig. 8(a) (c). Here, L = 10 8 and L = 10 3 are selected as ex- 
amples to compare the accuracy of approximation affected by the appearance 
of width one in rectangles. The discrepancy between the actual distribution 
and our approximation at t — 500 step in Fig. 9(a) (c) becomes larger than 
the corresponding results at t = 50 in Fig. 8(a) (c) for L = 10 3 , while the two 
distributions almost coincide in both Figs. 8 and 9 (b)(d) for L = 10 s . Figures 
8(e) (f) and 9(e) (f) show reasonable fittings of the actual distribution with the 
estimated lognormal functions in the cumulation for the body over the whole 
range of log 10 A and the tail restricted in the range of (e) log 10 A > 5 and (f ) 
Iog 10 A > 14. For the body and the tail, Figure 10 shows the lognormal den- 
sity functions. The short segment represents the slope of the exponent —2 for 
fragment-mass of glass [6] and areas enclosed by the city roads [8] [9] . The linear 
part in Fig. 10 seems to be longer for larger L further to the right, as similar to 
the area distributions of extreme rectangles in Fig. 2. Remember that a large 
L gives fine-grained divisions. This phenomenon may be consistent with the 
enhancement of power-law property as high impact in the fragmentation of 
glass [6], however the behavior in our model is not exactly the same because of 
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Fig. 8. (Color online) CCDs P(> A) at t = 50 step. The (red or green) solid and 
(black) dashed lines show our approximation by using (a)(b) the Poisson distribution 
of Eq.(10) or (c)(d) the solution of difference equation(5) and the averaged results 
by 100 samples of the actual divisions of rectangles, respectively. The (blue and 
cyan) lines show (e)(f) CCDs of the estimated lognormal functions for the body 
over the whole range of log 10 A and the tail in (e) log 10 A > 5 and (f) log 10 A > 14. 



the remaining lognormal property whose distribution rather resembles to the 
ones for road networks in German cities [10]. We note that the approximation 
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(e) lognormal, L = 10 3 (f) lognormal, L = 10 s 

Fig. 9. (Color online) CCDs P(> A) at t = 500 step. Approximation by the mix- 
ture YI1P1921 using (a)(b) Poisson distribution of Eq.(10) or (c)(d) the solution of 
differential Eq.(5) as pi, and (e)(f) the estimated lognormal functions for the body 
over the whole range of log 10 A and the tail in (e) log 10 A > 5 and (f) log 10 A > 14. 
The linetypes are the same as in Fig. 8. 

is no longer accurate for L = 10 8 in t 3> 500. The maximum difference < 0.1 
in Fig. 7 may be a criterion for whether it gives a good fitting or not. 
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(c) Tail, t = 50 (d) Tail, t = 500 

Fig. 10. (Color online) Approximation of area distributions estimated by lognormal 
functions in Fig. 8(e) (f) and Fig. 9(e) (f). The short segment guides the slope of -2. 

5 Conclusion 



Fat-tail distributions are pervasive in nature, and also appear in spatial net- 
works. In particular, lognormal and power-law distributions are familiar in 
fragments of glass [6,7], crack patterns, and areas enclosed by city roads 
[8,9,10]. Beyond the details of physical phenomena, it is useful to consider 
a common generation mechanism on the division process. Thus, we have in- 
vestigated a simple model for producing such distributions of areas enclosed by 
edges in a spatial network, which is an extension of MSQ networks [21,22,23] 
from the iterative divisions of squares to ones of rectangles. 

The stochastic division process makes a Markov chain in the random selec- 
tions of a rectangle face and of the division axes from the initial configuration 
of an L x L square. We have derived the exact solution of distribution at 
T m ax for the extreme rectangles with no more divisible edge(s) of width one 
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on a combinatorial analysis for a small L. It is also pointed out that, in the 
absorbing Markov chain, the iterative divisions of rectangles are equivalent to 
directional random walks with splitting into four particles. For a large L and 
t <C T max , we have discussed the distribution of areas on a continuous ap- 
proximation. As the distributions of layers and of areas restricted in a layer, 
we decompose the original distribution into two functions and consider the 
mixture of them, which corresponds to the Poisson and the gamma distribu- 
tions in Eqs. (10) and (11). In addition, we obtain the distribution of layers by 
the difference equation (5) more accurately. Simulation results show a good 
agreement of our approximation with the actual distribution in the divisions, 
and give a condition for the fitting. We obtain more accurate results, as the 
size L is larger and the time step t is smaller, since the layer of a face tends 
to be shallower and the effect of width one in rectangles becomes weaker. We 
emphasize that the decomposition into two distributions of layers and of areas 
restricted in a layer will be useful for investigating other phenomena, such as 
in broken fragments of glass and areas enclosed by city roads, with additional 
information of layers defined by the time sequence of divisions. 

We also confirm a slightly better agreement of our approximation for the 
m = 2 divisions which correspond to a crack model, though the meaning of 
network construction by bridgings should be discussed further from an appli- 
cation point of view. Unfortunately, relations to the mathematical properties 
of STIT tessellations [15] [16] [17] are unknown. In addition, there remain fur- 
ther challenges to the analysis for other properties, e.g., the lifetime of face 
or the distribution of edge lengths in our framework of stochastic processes, 
more rigorous investigation for the fitting functions (e.g., estimation by a dou- 
ble Pareto [3] [5]), extension to a preference selection model, and considering 
the division by any direction not limited to vertical and horizontal. 
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Appendix A: Fractal dimension 

We consider the fractal dimension df of a MSQ network [22] based on the self- 
similar tiling by squares at a finite time r < oo. Note that the asymptotical 
value is df — > 2 trivially for the infinite time r — > oo, since the division process 
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has no limitation and the whole two-dimensional space is embedded by the 
recursive subdivisions after a very long time. 

In general, for the number N b [l] of covered boxes by a measure e = 2~ l , it is 
defined as 

logiV fc [/] log N b [l] 

df = hm - — . . . = hm - — : — — . 
z^oolog(l/e) i^oo / -log 2 

For each visible level, the number N b [l] is counted as 

N b [2] = 2xN b [l}+n 2 (r), 
N b [3] = 2 x N b [2] +n 3 (r), 



where, in the above right-hand sides, the first term is due to the doubling 
of measure, and the second term comes from the one-to-one correspondence 
between four faces and cross edges generated by the quadratic division (e.g., 
in a clockwise mapping). The recursion is N b [l] = 2 l ~ 1 N b [l] + Y!i=2^ l ~ %n i{ T ) 
with N b [l] = 12 in an initial configuration of four squares. By substituting 
Eq.(8) with m = 4 into Yh=2 2 l ~ in i( T ), we derive 

yZ ol-iAi-l r'- 1 - T _ nl (2r) i ~ 1 - T 

^i=2 z ^ ~~ ^ i=2 2(i-l)! C ' 

= 2'(e 2r - l-Res)e- T /2, 

where Res denotes the residual for the higher-terms than I in the Taylor 
expansion of e x . 

For I ^> 1, we obtain 

\ogN b [l] ~ log(2 / (6 + e r /2)) ~ / • log 2 + r. (12) 

On the other hand, by using a generating function F t (z) == X^eLeauest z ^ w \ it 
is also represented as 

jV 6 [/]=4x ^ 2 Hw| = 4 x 2'F t (l/2), 

wdLeavest 

where | w | denotes the depth of face w, and Leaves t is a set of faces at the di- 
vision of i-step in the random quadtree, which corresponds to our hierarchical 
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network model. From the LEMMA 7.1 in [20], the expectation is 
£(JV t P]) = 4 x 2<£(F,(l/2)), 

W(z ))<exp ((4,-1)1^). 
When we set 3k + 1 = x, 

/3t-2 \ 

£(F t (l/2)) < exp ^ | — cfe) < (3t+ l) 1 / 3 . 

From the relation 1 + 3t = e 3r of the total number of faces, we obtain 

log£(N b [l}) <log4 + Z-log2 + -log(l + 3*) ~/-log2 + r. 

o 

This is equivalent to Eq.(12). 

We consider a typical time (r) = J2 T T Pi( r ) x m f° r the layer /, 

v ' 

where the factor m is due to the normalization, 

Thus, as an upper-bound for a large I, we obtain 

Z -log 2 + Z/m 1 

d f — -J— = 1 + — - = 1.36067. 

; Z-log2 m-log2 

Note that this value is between Koch curve: log 4/ log 3 = 1.26186 and Sierpin- 
ski carpet: log8/log3 = 1.89278 or Sierpinski gasket: log3/log2 = 1.58496. 

Appendix B: Extended preference model 

When a rectangle face is chosen (not uniformly at random but) proportionally 
to the power 7' of its depth / with a real parameter 7 > 0, we can approxi- 
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mately derive the distribution pi = ni/Af. Note that, in the selected division, 
a larger face (at a shallower layer) is preferred for 7 < 1, while a smaller face 
(at a deeper layer) is preferred for 7 > 1 [20]. 

The number rii of Z-th faces at time r follows 



_i = m i_ nj _ 1 _^ nj> /> 2 , (13) 

(14) 



As 7 ~ 1, we assume that 1/C is a constant. The reason is discussed later. By 
a variable transformation r = C^ l T, we rewrite Eq.(13) as 

7' drii drii •y 1 ( m \ 



Thus, we obtain the same form of Eq.(6) at 7 = 1 as follows 

dni m 

—— = — n;_i — rii. 
dT 7 

The solution is 



7; (/-i)! v-J c-i)! 



Here, we evaluate the assumption of a constant 1/C. Since the number of faces 
increases by m — 1 =3 (add four divisions, and delete one chosen face) at each 
time step t, the total number of faces is J2 rii = 1 + (m — l)t. It must be equal 
to Eq(9), so that = e^ m ~^ T . By applying dj^ni/dr = (m — l)Z)rij, the 
following relation must hold from Eqs.(13) and (14), 

1=1 ar 1=1 



therefore C = (J2il ln i) e ^ m ^ T ■ Indeed, around 7 ~ 1, 1/C is almost constant 
as approximated in Table 1. 
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Table 1 

The estimated values of 1/C by the Newton-Raphson method for varying a pa- 
rameter 7. Note that r takes a small value even for a huge network because r is 
a logarithmic timescale for a linear timescale t = 1,2, .. . of discrete steps in the 
relation 1 + (m - l)i = e (m " 1)r . 
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